Synchronization and Dephasing of Many-Body States in Optical Lattices 
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We introduce an approach to describe quantum-coherent evolution of a system of cold atoms in 
an optical lattice triggered by a change in superlattice potential. Using a time-dependent mean field 
description, we map the problem to a strong coupling limit of previously studied time-dependent 
BCS model. We compare the mean field dynamics to a simulation using light-cone methods and 
find reasonable agreement for numerically accessible times. The mean field model is integrable, and 
gives rise to a rich behavior, in particular to beats and recurrences in the order parameter, as well 
as singularities in the momentum distribution, directly measureable in cold atom experiment. 
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The approach to thermal equilibrium has been a cen- 
tral problem in statistical mechanics since Boltzmann's 
H-theorem. Recent advances in ultracold atoms make 
it possible to probe this approach in quantum-coherent 
many-body systems, due to the ability to change inter- 
actions in optical lattices [l], on a fast timescale. 

These new experimental opportunities stimulated the- 
oretical work on quantum dynamics in many-body sys- 
tems. Improvements in simulation algorithms, such as 
time-dependent density-matrix renormalization group [1, 
4] and light-cone methods [B| , make it possible to study 
these systems numerically for short times. However, the 
increase in entanglement entropy limits the simulation 
time [f|, calling for the development of new approaches. 

One simple-to-realize way to start a system out of equi- 
librium is to begin with an additional period-two modu- 
lation in a translationally invariant system. One case of 
this is an XXZ spin chain started from a Neel state at 
large Ising coupling [6i] . Another case proposed recently is 
a Bose gas in an optical lattice, with a period-two super- 
lattice initially superimposed[7|, causing the system to 
begin in a state with alternating filled and empty sites. 
The superlattice is then removed, and the system evolves 
under Bose-Hubbard dynamics. 

In this article we study interacting spinless fermions in 
a one-dimensional lattice, described by the Hamiltonian 

H = A ^2 (a]ai + i +h.c.J + ^ ^in i+ i, (1) 



i=l. .JV 



i=l...N 



where A is the hopping amplitude, hi — a\ai — ^. The ini- 
tial state, taken to be alternating filled and empty sites, 
hi = ±i, is created by an additional period- two potential 
that is removed at t = 0, after which the system evolves 
under H. The Hamiltonian (TJJ) describes the regime in 
which multiple occupancy of lattice sites is inhibited by 
repulsive interaction and/or the Pauli principle. 

To understand the evolution governed by |T]) we em- 
ploy a mean-field description of the problem (JTJ which 
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FIG. 1: Buildup of singular momentum distribution of pseu- 
dospin S%, a semicircle with the edge at k c — arcsinA. 
A series of traces computed for A = 0.5 are shown at times 
t = 0, 25, 50, 175, eventually converging to Eq. ([TTJ. Inset: 
the corresponding limiting fermion momentum distribution 
(n(k)} = {a{a k ) Q. 



uses the staggered density 
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as an order parameter. We show that the resulting mean- 
field dynamics is mathematically equivalent to time- 
dependent BCS dynamics @, @, 0, El, with, however, 
very different initial conditions. Comparison to the nu- 
merical results [j| for the Jordan- Wigner-equivalent XXZ 
spin chain is used to test validity of the mean-field ap- 
proach. 

We find that instead of simple relaxation to a steady 
state, the order parameter time evolution exhibits re- 
vivals. These revivals are understood as resulting from 
a buildup of singularities in the fermion momentum dis- 
tribution, illustrated in Fig[TJ The formation of a steady 
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state with a singular momentum distribution can be di- 
rectly tested in a free flight imaging experiment 13 1. 

The mean-field approximation can be constructed by 
replacing hih i+ i rs — 2 ( o Tr (t) J^^— l) l n i in the fermionic 
Hamiltonian |T]), which gives 
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(without loss of generality we set the hopping amplitude 
to A = |). In the evolution governed by H m { the quan- 
tity p v (t) is determined self-consistently via Eq.([2]). 

To establish equivalence of the problem ([3]) to the 
BCS dynamics [E S M, El we consider a translation 
invariant system on a ring, with momentum states |fc), 
— 7r < k < 7r. Initially, the system is half-filled and each 
pair of momentum states, |fc) and \k + n), contains ex- 
actly one particle. Momentum-nonconserving terms in 
the Hamiltonian ([3]) couple momenta k and k±it. Thus, 
at all times we can write the many-body wavefunction in 
a BCS-like form 

*(*)= II («*(*)<** + « fc (t)o t (fc))|0), (4) 

-7r/2<fe<7r/2 

where |0) is the vacuum state. The initial conditions are 
Uk (0) = Vk (0) = ;4= . The evolution equations are 

d t u k (t) = icos(k)u k (t) + 2i\pi r (t)v k (t), (5) 
d t v k (t) = -i cos(k)v k (t) + 2iXp 7T (t)u k (t), 

with the self-consistency condition |J2]) taking the form 
1 



N/2 



V He(B fc (t)t; fc (t)). 



(6) 



-7r/2<fc<ir/2 



Note that there are iV/2 different values of k in the sum 
in Eq. ©, so that -1/2 < p„(t) < 1/2. 

Given a pair, Uk,Vk, we define pseudo-spins [12j] by 
Su + iSu = UuVk ■ In terms 



of these classical variables the Hamiltonian reads 



S% = 



H s = - 



-7r/2<fe<7r/2 



2cos(fc)S*f 



w j:sisi, (?) 

' k,k' 



which, together with the usual Poisson brackets, 
{SajSb} = iSabcSc, reproduces the d yna mics ([5]). 

The canonical BCS Hamiltonian 12j differs from ([7]) 
in one important respect: the BCS problem has an ad- 
ditional coupling S^S^,. We circumvent this problem in 
two steps. First, we extend the Hamiltonian ([7]) to a 
twice larger momentum range — n < k < tt, doubling the 
number of A;-states. Simultaneously we replace the cou- 
pling 2A/ (N/2) by 2X/N. Consider any pair of momenta, 
k and k + tt; these two pseudospins feel opposite z fields 
and the same x field, so they evolve as 



Sk+Tr(t) — Sk\t), 



stUt) 



-st' z (t), 




Time 



Time 



FIG. 2: Comparison of mean-field (black, Eq.©) and exact 
dynamics (red, Eq.©) simulated using light-cone methods [jj. 
The revival in a) at t > 20 can be modeled by Eq. pOjl , 



where individual spin polarizations for —tt/2 < k < n/2 
evolve in the same way as in the original problem (UJ . 

After the states are doubled, the net y polarization 
J2 k taken over the extended range — tt < k < tt van- 
ishes at all times, while the net x polarization S k 
does not change, and so we can then add the interaction 
S^Sf, to the original Hamiltonian ((JJ without changing 
the dynamics. Thus, we arrive at the BCS Hamiltonian: 
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— 7T<k<7T 
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2coskS z k + -J2 S kS X k'+S y k Sl. (8) 

k.k' 



With the initial state S% = 1, Sf z = 0, the problem © 
yields the dynamics of p-„(t) identical to (JTJ) . 

As a sanity check of the mean-field approximation we 
use comparison to the simulation of the XXZ spin- i chain 
with the zz coupling of strength A: Hxxz — Y^,i S?Sf +1 + 
SfSf +1 + \SfS* +1 which is Jordan- Wigner-equivalent to 
((T|). The initial conditions are taken to be the Neel state: 
alternating spin up and down. The staggered density © 
translates to the Neel order parameter 



M(t) 
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In Fig. [2ja) we compare the XXZ simulation done us- 
ing the light-cone method [f| to the mean-field dynamics 
(IS]) for A = 0.5. The mean-field dynamics was simulated 
using 4th order Runge-Kutta with 25000 pairs of modes 
Uk, Vk with < k < n/2 and a time step of 0.025; the 
parameters were chosen to assure insensitivity of the re- 
sults to finite size and finite timestep. The behavior of 
the exact M(t) and the mean- field p-n(t) is quite simi- 
lar, with M(t) oscillating slightly faster than p^(t). The 
later times of the light-cone dynamics show noise from 
sampling errors. Note that at times of just past 20, the 
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mean-field dynamics shows a revival: the amplitude of 
oscillations begins to increase again. This revival is not 
far outside the times reached with the light-cone meth- 
ods, and may be accessible with more numerical effort. 

As we will see below, for A > f the mean-field dynam- 
ics predicts a non- vanishing asymptotic value p n (t — ► oo), 
while simulation of the XXZ chain indicates that p n {t) 
rapidly decays to zero. Nevertheless, as illustrated in 
Fig. [2]Jb) , even at A = 1 there is a reasonably good agree- 
ment between the two approaches at short times, with 
both M(t) and p n (t) decaying faster than for A = 0.5. 

Turning to discuss different regimes, we note that be- 
cause both the Hamiltonian ([3]) and the initial values 
Ufc(0) and Wfc(0) are real, the mean-field dynamics of Pn-(t) 
is time reversal invariant. Further, the initial state is in- 
variant under the orthogonal operator O = J} fc S¥, which 
anti-commutes with the first term in H m { and commutes 
with the second. Combining these two statements, we see 
that the dynamics depends only on the magnitude of A 
and not on its sign (all this is also true for the XXZ chain 
started in the Neel state, where we set O = Hiodd ^i)- 

In the absence of interaction, A = Ojthe decay of p-„(t) 
follows cos(2i + 5)/t 1 / 2 with 8= f [iQ. At < |A| < 1, 
the system is in a "dephasing" regime, and p w is well- 
described at long times by a sum of two frequencies beat- 
ing together with a power-law envelope decaying as t~ 3 / 2 : 

p„(t) ~ (aiCos(wii + 5i) + a2Cos(w2£ + £2))/i 3/2 - (10) 

We will see that this unexpected behavior, leading to 
revivals in /?7r(i), signals formation of a singularity in 
k space of the asymptotic polarization SI. Below, we 
use integrals of motion of the BCS dynamics to show, in 
agreement with numerics, that uji = 2, u>2 = 1\J\ — A 2 . 

One can expect that strong interaction will stabilize 
the state with density modulation ([2]) . In agreement with 
this intuition, we find that for |A| > 1 the system is 
in the "polarized" regime, with non- vanishing p^(t — > 
oo) = (1/2)^1 — 1/A 2 . We obtain this asymptotic value 
analytically, and confirm it numerically. The approach to 
the asymptotic value is described by p n (t) — ^/l — 1/A 2 ~ 
cos(wt + J)/^ 1 / 2 . At the phase transition at A = 1, we 
observe numerically that p^it) ~ cos(2< + S)/t 3 ^ 2 . 

The behavior in the dephasing regime can be qualita- 
tively understood as follows. In the non-interacting case, 
all of the pseudo-spins remain in the x — y plane, with 
S%(t) = 0. Initially, they all point in the x-direction, 
and dephase over time, leading to the 1/yi decay in p^. 
When interaction is turned on, the spins begin to move 
out of the x — y plane and polarize in the z-direction. The 
appearance of a net z-polarization is required by conser- 
vation of energy: as the ferromagnetic S X S X term in the 
energy §S§ is becoming less negative because of dephas- 
ing, the S z contribution must become more negative. 

However, the asymptotic distribution of S z is not the 
thermal distribution. Instead, it shows a square-root sin- 



gularity at k = k c = arcsinA (see Fig. [IJ. We find, ana- 
lytically and numerically (see Figfl]), that Sjj.(t — > oo) = 
for |fc| > k c , whereas 



d) = i\/cos 2 A: -1 + A 2 , |fc|<fc c . (11) 



The van Hove singularity at the band edge and the sin- 
gularity at k c give rise to the frequencies u>i, ll>2 in (| 10[) . 

An analytic insight into the behavior in the dephasing 
regime can be gained as follows. The BCS dynamics has 
infinitely many commuting integrals of motion that can 
be written as an energy dependent Lax vector Q 



L(£) = z + 2A^ 



Sdk 



where in our case £ = cosk and ...S^/ = J 

The asymptotic polarization can be found by compar- 
ing the values L(£) in the initial and asymptotic states 
The initial state, polarized in the x direction, gives 



L 2 (0 = L 2 1+ Ll 



l + A 2 /(£ 2 -l). 



(12) 



In the asymptotic state, because the x and y components 
are dephased, we can approximate: 



L 2 (0^L 2 = (l + 2A ^ 



S z k dk 



(13) 



(£ — cos/c)27r / 

Comparing (|12[) and Q13p we obtain an integral equation 



2A 



S%dk 



-ir (£ — COs/c)27T 



/g 2 -l + A 2 



(14) 



where ^ is treated as a complex variable, Im£; > 0. 
Changing variable to x = cos k, we rewrite (fl4l) as 



J_ f +1 f(x)dx 
2tt 



-l 



X 



l£ 2 ~ 1 + A 2 



(15) 



where f(x) = AXS z {x)/VT — x 2 is the unknown function. 
This equation can be solved using Cauchy's formula, by 
writing f(x) — f+{x) the functions f±{x) being 

analytic in the upper/lower complex halfplane, respec- 



tively. The result is f(x) = 2 Im y ^-prjp- , which yields 
the semicircle dependence (fTTj) . 

The t -3 / 2 power law envelope in (fT0|) is more difficult 
to understand. Consider the band edge singularity at 
k = 0. At A = 0, this gives rise to a contribu- 
tion to due to a coherent contribution of the modes 
with k < \j\ft. For A > 0, the x — y component of the 
spins near k — is proportional to k, which suggests in a 
scaling picture that S X (Q) should decay with an envelope 
of l/Vt, which we do observe numerically. Thus, if the 
order parameter were due to a coherent contribution of 
spins with k < we would expect a. 1/t decay in p„, 
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FIG. 3: The quantity S? um (fc) = £ E < fc << fc Sfr. 
(black) and A = 0.5 (red). Snapshots at t — 400 are shown, 
with y-axis multiplied by an arbitrary scaling factor in both 
cases. 

rather than 1/t 3 / 2 as observed. To understand this bet- 
ter, we analyze the cumulative order parameter S^ um (k), 
plotting it in Fig. [3] for a particular snapshot in time. It 
can be seen that the contribution of the low-lying modes 
is largely canceled out by higher modes in the interacting 
case, while the contribution of the higher modes averages 
out in the non-interacting case. 

To identify the boundary of the dephasing regime, we 
recall that in the BCS problem the \arge-t asymptotic is 
governed by the complex roots of the spectral equation 
L 2 (£) = d, [l(J In our case, for the initial state 
polarized in the x-direction, the spectral equation can be 
factorized as L 2 = (L3 + iLi)(L^ — iL\) = 0, giving 

- A • (16) 



1 ± i2\ 



4-7r(£ — cos k) 



= 1 ± 



= 



This equation has complex imaginary roots £ = 
±iy/l — A -2 if A > 1, and has no complex solutions for 
A < 1. Thus, for A > 1, we have a non-zero asymp- 

~ 2 which we confirmed by 



totic value of — , 
numerical simulations. 

We have tried non-integrable deformations of the 
model, by making the S X S X coupling between modes 
weakly dependent on k. Even a very small change re- 
moves the singularity at k c , though a smooth kink re- 
mains, causing the contribution to p v {t) at frequency 
L02 — 2 cos k c to decay exponentially in time, rather than 
as a power law. At early times, beats are still observable. 

Now we briefly discuss application of our results to 
bosons in an optical lattice, described by [HI 



(17) 



As above, the initial state of alternating filled and empty 
sites is imposed by an additional period-two potential. 



Because the mean-field dynamics with the staggered 
density order parameter p n (t) — ~ 1/2) 

does not depend on statistics, the bosonic and fcrmionic 
Hamiltonians (jTTJ) , ([1]) yield the same mean-field evolu- 
tion. In fact, in the bosonic case, for any initial con- 
ditions where we alternate a site with n particles with 
an empty site, we obtain the same mean-field equations 
up to rescaling U by dividing it by n. We expect the 
mean-field to become more accurate for larger n. 

As a result, the mean- field theory for the bosonic 
system predicts the asymptotic momentum distribution 
n(k) — (btbf.) which is similar to fermionic n(k) (see 
Fig. [I]). This momentum distribution can be measured in 
a cold atom experiment using time-of-flight. 

Unlike the fermion problem (p}, the Bose-Hubbard 
model (jTTJ) is not integrable, and therefore at long times 
the momentum distribution of the particles should ther- 
malize. Still, at short times the kink in the distribution 
n(k) may be observable. If present, it will manifest itself 
also in revivals of the staggered density amplitude p„(t). 

As the interaction U increases, it is no longer valid 
to use mean-field theory. However, because at U = 00 
the bosonic problem reduces to non- interacting fermions, 
for large U we can use second order perturbation theory 
to map the problem onto a system of hard-core bosons 
with weak attractive interactions ^- Ei + 

TrLi b\bib\ +x bi-x + h.c. For large U, both terms are 
now weak and can be treated by mean-field theory. The 
treatment of the first term is as before, whereas the sec- 
ond term leads to a momentum dependent coupling in 
the BCS mean-field which breaks integrability. 

In contrast to that, the fermion problem is integrable, 
and thus its dynamics is not be ergodic. Thus our main 
results, the revivals in the order parameter and the for- 
mation of a singular momentum distribution, may persist 
in the fermion case even at the times longer than those 
described by our mean-field approach. 
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